Load packages

## add 'developer' to packages to be installed from github
x <- c("remotes", "lubridate", "readxl", "pbapply", "viridis", "ggplot2", "kableExtra", "knitr", "formatR", "randomForest", "Rtsne", "MASS", "adehabitatHR", "sp", "GGally", "spatstat", "raster", "vegan", "maRce10/PhenotypeSpace")

aa <- lapply(x, function(y) {
  
  # get pakage name
  pkg <- strsplit(y, "/")[[1]]
  pkg <- pkg[length(pkg)]
  
  # check if installed, if not then install 
  if (!pkg %in% installed.packages()[,"Package"])  {

      if (grepl("/", y))  remotes::install_github(y, force = TRUE) else
    install.packages(y) 
    }

  # load package
  a <- try(require(pkg, character.only = T), silent = T)

  if (!a) remove.packages(pkg)
  })

Functions and global parameters

opts_knit$set(root.dir = "..")

# set evaluation false
opts_chunk$set(fig.width = 6, fig.height = 3, warning = FALSE, message = FALSE, tidy = TRUE)

read_excel_df <- function(...) data.frame(read_excel(...))

# for reading months in english format
sl <- Sys.setlocale(locale = "en_US.UTF-8")
# read ext sel tab calls
sels <- read.csv("./data/processed/tailored_budgie_calls_sel_tab.csv")

# keep only spectrographic parameters
sels <- sels[, c("sound.files", "selec", "duration", "meanfreq", "sd", "freq.median",
    "freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom",
    "mindom", "maxdom", "dfrange", "modindx", "meanpeakf")]

sels$ID <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][1])
sels$month <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][2])
sels$day <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][3])
sels$year <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][4])

sels$date <- paste(sels$day, substr(sels$month, 0, 3), sels$year, sep = "-")

sels$date <- as.Date(sels$date, format = "%d-%b-%Y")

Counts per individual

df <- as.data.frame(table(sels$ID))

names(df) <- c("ID", "Sample_size")

df <- df[order(df$Sample_size, decreasing = FALSE), ]

kb <- kable(df, row.names = FALSE)

kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))

print(kb)
ID Sample_size
132YGMM 6
125YGMM 12
160YGHM 16
310YGHM 21
393YGHM 25
303YGHM 37
398WBHM 41
365YGHM 44
399YGLM 46
300YGHM 47
270YGHM 64
407YGHM 97
386WBHM 100
377WWLM 108
367WWNM 119
371YYLM 149
397YGHM 175
378YYLM 195
404WBHM 196
258YGHM 223
389WWLM 230
262YGHM 306
400YGHM 360
388YGLM 377
327YYLM 404
396YBHM 455
403WBHM 566
356WBHM 761
361YGHM 767
402YGHM 849
395WBHM 854
360YGHM 900
390WBHM 935
405YBHM 975
385YBHM 1256
394YBHM 1693

Add metadata

metadat <- read_excel_df("./data/raw/INBREStress_MasterDataSheet_14Nov19.xlsx")

# head(metadat)

sels$ID[sels$ID == "125YGMM"] <- "125YGHM"
sels$ID[sels$ID == "394YBHM"] <- "394WBHM"

# setdiff(sels$ID, metadat$Bird.ID) setdiff(metadat$Bird.ID, sels$ID)

sels$treatment <- sapply(1:nrow(sels), function(x) {

    metadat$Treatment[metadat$Bird.ID == sels$ID[x]][1]

})


sels$treatment.room <- sapply(1:nrow(sels), function(x) {

    metadat$Treatment.Room[metadat$Bird.ID == sels$ID[x]][1]

})


sels$round <- sapply(1:nrow(sels), function(x) {

    metadat$Round[metadat$Bird.ID == sels$ID[x]][1]

})

sels$source.room <- sapply(1:nrow(sels), function(x) {

    metadat$Source.Room[metadat$Bird.ID == sels$ID[x]][1]

})

sels$record.group <- sapply(1:nrow(sels), function(x) {

    metadat$Record.Group[metadat$Bird.ID == sels$ID[x]][1]

})


# add week
out <- lapply(unique(sels$round), function(x) {

    Y <- sels[sels$round == x, ]

    min_date <- min(Y$date)
    week_limits <- min_date + seq(0, 100, by = 7)

    Y$week <- NA
    for (i in 2:length(week_limits)) Y$week[Y$date >= week_limits[i - 1] & Y$date <
        week_limits[i]] <- i - 1

    return(Y)
})

sels <- do.call(rbind, out)


sels$cort.baseline <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


sels$cort.stress <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


sels$weight <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


sels$breath.count <- sapply(1:nrow(sels), function(x) {

    if (sels$week[x] == 1)
        out <- metadat$D3.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 2)
        out <- metadat$D7.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 3)
        out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 4)
        out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    if (sels$week[x] == 5)
        out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]

    return(out)
})


# remove week 5
sels <- sels[sels$week != 5, ]

Exploratory graph

Physiological parameters

breath.count <- stack(metadat[, c("D3.Breath.Count", "D7.Breath.Count", "D14.Breath.Count",
    "D21.Breath.Count", "D28.Breath.Count")])

weight <- stack(metadat[, c("D3.Bird.Weight..g.", "D7.Bird.Weight..g.", "D14.Bird.Weight..g.",
    "D21.Bird.Weight..g.", "D28.Bird.Weight..g.")])

cort.stress <- stack(metadat[, c("D3.CORT.Stress", "D7.CORT.Stress", "D14.CORT.Stress",
    "D21.CORT.Stress", "D28.CORT.Stress")])

cort.baseline <- stack(metadat[, c("D3.CORT.Baseline", "D7.CORT.Baseline", "D14.CORT.Baseline",
    "D21.CORT.Baseline", "D28.CORT.Baseline")])



stress <- data.frame(metadat[, c("Bird.ID", "Treatment", "Round", "Treatment.Room")],
    week = breath.count$ind, breath.count = breath.count$values, weight = weight$values,
    cort.stress = cort.stress$values, cort.baseline = cort.baseline$values)

# head(stress)

stress$week <- factor(sapply(strsplit(as.character(stress$week), "\\."), "[[", 1),
    levels = c("D3", "D7", "D14", "D21", "D28"))

stress$day <- as.numeric(gsub("D", "", as.character(stress$week)))
stress$round <- paste("Round", stress$Round)

stress$treatment <- factor(stress$Treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

# remove 5th week
stress <- stress[stress$week != "D28", ]

stress_l <- lapply(stress$Bird.ID, function(x) {
    X <- stress[stress$Bird.ID == x, ]

    X$breath.count <- X$breath.count - X$breath.count[X$week == "D3"]
    X$weight <- X$weight - X$weight[X$week == "D3"]
    X$cort.stress <- X$cort.stress - X$cort.stress[X$week == "D3"]
    X$cort.baseline <- X$cort.baseline - X$cort.baseline[X$week == "D3"]

    return(X)
})

stress <- do.call(rbind, stress_l)

ggplot(data = stress, aes(x = day, y = weight, group = Bird.ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Weight (g)", x = "Day sample was taken") +
    theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
    labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

ggplot(data = stress, aes(x = day, y = breath.count, group = Bird.ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Breath count", x = "Day sample was taken") +
    theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
    labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

ggplot(data = stress, aes(x = day, y = cort.baseline, group = Bird.ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Baseline CORT", x = "Day sample was taken") +
    theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
    labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

ggplot(data = stress, aes(x = day, y = cort.stress, group = Bird.ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Stress CORT", x = "Day sample was taken") +
    theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
    labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

Behavioral parameters

Call counts per treatment

call_count <- aggregate(selec ~ week + treatment + round, data = sels, length)
call_count$Week <- factor(call_count$week, levels = (unique(call_count$week))[5:1])

call_count$round <- paste("Round", call_count$round)

call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

ggplot(data = call_count, aes(x = treatment, y = selec, fill = Week)) + geom_bar(stat = "identity") +
    scale_fill_viridis_d(begin = 0.1, end = 0.9) + coord_flip() + facet_wrap(~round,
    ncol = 2) + labs(y = "Number of calls", x = "Treatment") + theme_classic(base_size = 30) +
    theme(legend.position = c(0.8, 0.15))

Call counts per bird and treatment

# Call counts per treatment by individual (in color)

call_count <- aggregate(selec ~ week + treatment + round + ID, data = sels, length)
call_count$Week <- factor(call_count$week, levels = (1:5))

call_count$round <- paste("Round", call_count$round)

call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))


ggplot(data = call_count, aes(x = week, y = selec, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Stress CORT", x = "Day sample was taken") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Acoustic parameters

call_week_params <- aggregate(cbind(duration, meanfreq, sd, freq.median, freq.IQR,
    time.IQR, skew, kurt, sp.ent, time.ent, entropy, meandom, mindom, maxdom, dfrange,
    modindx, meanpeakf) ~ week + treatment + round, data = sels, mean)

call_week_params_sd <- aggregate(cbind(duration, meanfreq, sd, freq.median, freq.IQR,
    time.IQR, skew, kurt, sp.ent, time.ent, entropy, meandom, mindom, maxdom, dfrange,
    modindx, meanpeakf) ~ week + treatment + round, data = sels, sd)

names(call_week_params_sd) <- names(call_week_params) <- c("week", "treatment", "round",
    "Duration", "Mean_frequency_(kHz)", "SD_frequency_(kHz)", "Median_frequency_(kHz)",
    "Interquantile_frequency_range_(kHz)", "Interquantile_time_range_(s)", "Skewness",
    "Kurtosis", "Spectral_entropy", "Time_entropy", "Overall entropy", "Mean_dominant_frequency_(kHz)",
    "Minimum_dominant_frequency_(kHz)", "Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)",
    "Modulation_index", "Mean_peak_frequency_(kHz)")


call_week_params_sd <- call_week_params_sd[, c("Duration", "Mean_frequency_(kHz)",
    "SD_frequency_(kHz)", "Median_frequency_(kHz)", "Interquantile_frequency_range_(kHz)",
    "Interquantile_time_range_(s)", "Skewness", "Kurtosis", "Spectral_entropy", "Time_entropy",
    "Overall entropy", "Mean_dominant_frequency_(kHz)", "Minimum_dominant_frequency_(kHz)",
    "Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)", "Modulation_index",
    "Mean_peak_frequency_(kHz)")]

names(call_week_params_sd) <- paste(names(call_week_params_sd), "sd", sep = ".")


call_week_params <- cbind(call_week_params, call_week_params_sd)
call_week_params$round <- paste("Round", call_week_params$round)
call_week_params$treatment <- factor(call_week_params$treatment, levels = c("Control",
    "Medium Stress", "High Stress"))


pd <- position_dodge(0.3)

ggs <- lapply(c("Duration", "Mean_frequency_(kHz)", "SD_frequency_(kHz)", "Median_frequency_(kHz)",
    "Interquantile_frequency_range_(kHz)", "Interquantile_time_range_(s)", "Skewness",
    "Kurtosis", "Spectral_entropy", "Time_entropy", "Overall entropy", "Mean_dominant_frequency_(kHz)",
    "Minimum_dominant_frequency_(kHz)", "Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)",
    "Modulation_index", "Mean_peak_frequency_(kHz)"), function(i) {
    call_week_params$var <- call_week_params[, i]
    call_week_params$var.min <- call_week_params$var - call_week_params[, paste(i,
        "sd", sep = ".")]
    call_week_params$var.max <- call_week_params$var + call_week_params[, paste(i,
        "sd", sep = ".")]

    ggplot(call_week_params, aes(x = week, y = var, col = treatment)) + geom_point(size = 4,
        position = pd) + geom_errorbar(aes(ymin = var.min, ymax = var.max), width = 0.3,
        position = pd, size = 1.7) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
        facet_wrap(~round, ncol = 2) + labs(y = i, x = "Week") + theme_classic(base_size = 30) +
        theme(legend.position = c(0.8, 0.25))
})

ggs
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

Acoustic space projection

scale_acous_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median",
    "freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom",
    "mindom", "maxdom", "dfrange", "modindx", "meanpeakf")])

urfmod <- randomForest(x = scale_acous_param, importance = TRUE, proximity = TRUE,
    ntree = 10000)

saveRDS(urfmod, "./data/processed/unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")

mds_urf_prox <- cmdscale(1 - urfmod$proximity, k = 2)

saveRDS(mds_urf_prox, "./data/processed/MDS_unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")

Random forest

urfmod <- readRDS("./data/processed/unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")

mds_urf_prox <- readRDS("./data/processed/MDS_unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")

sels$MDS1 <- mds_urf_prox[, 1]
sels$MDS2 <- mds_urf_prox[, 2]

# print(1- sum(urfmod$votes[,1] > urfmod$votes[,2])/ nrow(urfmod$votes))

Proportion of real data classified as synthetic data by the unsupervised random forest: 0.00506

ggplot(sels, aes(x = MDS1, y = MDS2, col = as.factor(round))) + geom_point() + labs(color = "Round") +
    scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 30) + theme(legend.position = c(0.62,
    0.3)) + guides(colour = guide_legend(override.aes = list(size = 10)))

vimp <- as.data.frame(vimp)
vimp$params <- rownames(vimp)
vimp <- vimp[order(vimp$MeanDecreaseAccuracy), ]
vimp$params <- as.factor(vimp$params)
vimp$params <- factor(vimp$params, levels = vimp$params[!duplicated(vimp$params)])
vimp <- vimp[(nrow(vimp) - 15):nrow(vimp), ]

ggplot(vimp, aes(x = MeanDecreaseAccuracy, y = params)) + geom_segment(aes(yend = params),
    xend = 0, color = "grey50", size = 1) + geom_point(size = 10, col = viridis(10,
    alpha = 0.6)[2]) + labs(x = "Mean decrease accuracy (RF)", y = "Predictors") +
    theme(legend.key.size = unit(2, "lines")) + scale_color_discrete(name = "Parameter set") +
    scale_fill_discrete(guide = "none") + theme_classic(base_size = 25)

PCA

pca <- prcomp(x = sels[, c("duration", "meanfreq", "sd", "freq.median", "freq.IQR",
    "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom",
    "maxdom", "dfrange", "modindx", "meanpeakf")], scale. = TRUE)

Y <- as.data.frame(pca$x[, 1:2])
names(Y) <- c("PC1", "PC2")

sels <- data.frame(sels, Y)

ggplot(sels, aes(x = PC1, y = PC2, col = as.factor(round))) + geom_point() + labs(color = "Round") +
    scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) + theme(legend.position = c(0.9,
    0.8)) + guides(colour = guide_legend(override.aes = list(size = 10)))

t-SNE

scale_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median", "freq.IQR",
    "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom",
    "maxdom", "dfrange", "modindx", "meanpeakf")])

tsne <- Rtsne(scale_param, dims = 2, perplexity = 30, verbose = FALSE, max_iter = 5000)

saveRDS(tsne, "./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
tsne <- readRDS("./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")

Y <- as.data.frame(tsne$Y)
names(Y) <- c("TSNE1", "TSNE2")

sels <- data.frame(sels, Y)

ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(round))) + geom_point() +
    labs(color = "Round") + scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) +
    theme(legend.position = c(0.955, 0.25)) + guides(colour = guide_legend(override.aes = list(size = 10)))

Acoustic space by individual (each point) for 3 the dimension reduction methods

  • Only individuals with at least 20 calls
  • Using rarefaction with n = the minimum sample size across all individuals
tab <- table(sels$ID)

Y <- sels[sels$ID %in% names(tab[tab > 20]), ]

pca_areas <- rarefact_space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID",
    parallel = 4, pb = TRUE, type = "mcp")

rf_mds_areas <- rarefact_space_size(X = Y, dimensions = c("MDS1", "MDS2"), group = "ID",
    parallel = 4, pb = TRUE, type = "mcp")

tsne_areas <- rarefact_space_size(X = Y, dimensions = c("TSNE1", "TSNE2"), group = "ID",
    parallel = 4, pb = TRUE, type = "mcp")

areas <- data.frame(pca_areas[, -ncol(pca_areas)], pca.area = pca_areas$mean.size,
    rf.mds.area = rf_mds_areas$mean.size, tsne.area = tsne_areas$mean.size)

saveRDS(areas, "./data/processed/acoustic_space_area_by_individual_for_each_dimension_reduction_method.RDS")
areas <- readRDS("./data/processed/acoustic_space_area_by_individual_for_each_dimension_reduction_method.RDS")

ggpairs(areas, columns = which(names(areas) %in% c("pca.area", "rf.mds.area", "tsne.area")),
    columnLabels = c("PCA", "Rand.for. MDS", "t-SNE")) + theme_minimal(base_size = 30)

Individual acoustic space across weeks

  • PCA used for dimensionality reduction
  • Only weeks with at least 6 calls were included
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")

tab <- table(sels$ID.week)

# sort(tab)

Y <- sels[sels$ID.week %in% names(tab)[tab >= 6], ]

# run with rarefaction
areas_by_week_raref <- rarefact_space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID.week",
    parallel = 4, pb = TRUE, type = "mcp")

# run without rarefaction
areas_by_week <- space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID.week",
    cl = 4, pb = TRUE, type = "mcp")

areas_by_week$raref.area <- areas_by_week_raref$mean.size

# add metadata
areas_by_week$ID <- sapply(strsplit(areas_by_week$group, "-"), "[[", 1)

areas_by_week$week <- factor(sapply(strsplit(areas_by_week$group, "-"), "[[", 2),
    levels = 1:5)

areas_by_week$treatment <- sapply(1:nrow(areas_by_week), function(x) metadat$Treatment[metadat$Bird.ID ==
    areas_by_week$ID[x]][1])

areas_by_week$round <- sapply(1:nrow(areas_by_week), function(x) metadat$Round[metadat$Bird.ID ==
    areas_by_week$ID[x]][1])

areas_by_week$round <- paste("Round", areas_by_week$round)

call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

areas_by_week$treatment <- factor(areas_by_week$treatment, levels = c("Control",
    "Medium Stress", "High Stress"))

saveRDS(areas_by_week, "./data/processed/acoustic_space_area_by_individual_and_week.RDS")

With rarefaction

areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")

ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

With rarefaction and compared to week 1

areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")

areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {

    X <- areas_by_week[areas_by_week$ID == x, ]
    X$raref.area <- X$raref.area - X$raref.area[X$week == min(as.numeric(as.character(X$week)))]
    return(X)
}))

ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Without rarefaction

ggplot(data = areas_by_week[areas_by_week$week != 5, ], aes(x = week, y = size, group = ID,
    col = treatment)) + geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1,
    end = 0.9) + facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Without rarefaction and compared to week 1

areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {

    X <- areas_by_week[areas_by_week$ID == x, ]
    X$size <- X$size - X$size[X$week == min(as.numeric(as.character(X$week)))]
    return(X)
}))

ggplot(data = areas_by_week, aes(x = week, y = size, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Acoustic space area by number of calls per individual and week

ggplot(data = areas_by_week, aes(x = n, y = size, col = treatment)) + geom_point() +
    # geom_smooth(size = 1.2, method = 'lm') +
scale_color_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~round, ncol = 2, scales = "free_x") +
    labs(y = "Acoustic space area", x = "Number of calls (n)") + theme_classic(base_size = 30) +
    theme(legend.position = c(0.8, 0.25))

  • Acoustic space area by number of calls per individual and week by treatment
  • Must have at least 2 calls
ggplot(data = areas_by_week, aes(x = n, y = log(size + 10), col = treatment)) + geom_point() +
    geom_smooth(size = 1.2, method = "lm") + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    # facet_wrap(~ round, ncol = 2, scales = 'free_x') +
labs(y = "log(acoustic space area)", x = "Number of calls (n)") + theme_classic(base_size = 30) +
    theme(legend.position = c(0.8, 0.7))

Acoustic space overlap

  • 2 methods: pairwise distance between observations and kernel density overlap
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 2], ]

# run with density
dens_overlap_by_id <- space_similarity(X = Y, dimensions = c("PC1", "PC2"), group = "ID",
    parallel = 10, pb = TRUE, outliers = 0.95, type = "mean.density.overlap")

dist_overlap_by_id <- space_similarity(X = Y, dimensions = c("PC1", "PC2"), group = "ID",
    parallel = 10, pb = TRUE, outliers = 0.95, type = "distance")

saveRDS(list(dens_overlap_by_id = dens_overlap_by_id, dist_overlap_by_id = dist_overlap_by_id),
    "./data/processed/acoustic_space_overlap_by_individual_2_metrics.RDS")

Mantel test between the 2 methods

attach(readRDS("./data/processed/acoustic_space_overlap_by_individual_2_metrics.RDS"))

mantel.rtest(as.dist(1 - rectangular_to_triangular(dens_overlap_by_id, distance = FALSE)),
    as.dist(rectangular_to_triangular(dist_overlap_by_id)), nrepet = 10000)
## Monte-Carlo test
## Call: mantelnoneuclid(m1 = m1, m2 = m2, nrepet = nrepet)
## 
## Observation: 0.6632477 
## 
## Based on 10000 replicates
## Simulated p-value: 9.999e-05 
## Alternative hypothesis: greater 
## 
##       Std.Obs   Expectation      Variance 
##  5.5337964986 -0.0001518025  0.0143715608

Individual overlap to its own group

sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 2], ]

group_ovlp_l <- lapply(unique(Y$ID.week), function(x) {

    # group data
    X <- Y[Y$record.group == Y$record.group[Y$ID.week == x][1] & Y$week == Y$week[Y$ID.week ==
        x][1], ]

    # label all other individuals as group
    X$ID2 <- ifelse(X$ID.week == x, "focal", "group")

    # run with density
    ovlp <- if (min(table(X$ID2)) > 4)
        space_similarity(X = X, dimensions = c("PC1", "PC2"), group = "ID2", parallel = 10,
            pb = FALSE, outliers = 0.95, type = "mean.mcp.overlap") else matrix(rep(NA, 3), nrow = 1)

    dsts <- if (min(table(X$ID2)) > 1)
        space_similarity(X = X, dimensions = c("PC1", "PC2"), group = "ID2", parallel = 10,
            pb = FALSE, outliers = 0.95, type = "distance") else matrix(rep(NA, 3), nrow = 1)

    out <- data.frame(ID = Y$ID[Y$ID.week == x][1], group = Y$record.group[Y$ID.week ==
        x][1], week = Y$week[Y$ID.week == x][1], overlap.to.group = ovlp[1, 3], distance.to.group = dsts[1,
        3], treatment = Y$treatment[Y$ID.week == x][1], round = Y$round[Y$ID.week ==
        x][1])

    return(out)

})

group_ovlp <- do.call(rbind, group_ovlp_l)

saveRDS(group_ovlp, "./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")
group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")


group_ovlp <- do.call(rbind, lapply(unique(group_ovlp$ID), function(x) {

    X <- group_ovlp[group_ovlp$ID == x, ]
    X$overlap.to.group <- X$overlap.to.group - X$overlap.to.group[X$week == min(as.numeric(as.character(X$week)))]
    X$distance.to.group <- X$distance.to.group - X$distance.to.group[X$week == min(as.numeric(as.character(X$week)))]
    return(X)
}))


group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to group acoustic space",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual distance to group acoustic space",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

# Partial mantel test between acoustic space
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")

tab <- table(sels$ID.week)

Y <- sels[sels$ID.week %in% names(tab)[tab >= 2], ]

treatment_labs <- sapply(unique(c(dens_overlap_by_id$group.1, dens_overlap_by_id$group.2)),
    function(x) Y$treatment[Y$ID == x][1])

treatment_mat <- binary_matrix(labels = treatment_labs)

round_labs <- sapply(unique(c(dens_overlap_by_id$group.1, dens_overlap_by_id$group.2)),
    function(x) Y$round[Y$ID == x][1])

round_mat <- binary_matrix(labels = round_labs)

# mantel.rtest(as.dist(treatment_mat), as.dist(1 - dist_overlap_by_id))

# mantel.partial(as.dist(treatment_mat), as.dist(1 - dens_overlap_by_id), zdis
# = as.dist(round_mat))
sels$ID.week <- paste(sels$ID, sels$week, sep = "-")

tab <- table(sels$ID.week)

Y <- sels[sels$ID.week %in% names(tab)[tab >= 5], ]

# run with rarefaction
overlap_by_week <- space_similarity(X = Y, dimensions = c("PC1", "PC2"), group = "ID.week",
    parallel = 1, pb = TRUE, outliers = 0.95, type = "mean.mcp.overlap")

saveRDS(overlap_by_week, "./data/processed/acoustic_space_overlap_by_individual_and_week.RDS")
trait_space_plot <- function(X, dimensions, group, indices, basecex = 1, title = "",
    colors = c("#3E4A89FF", "#35B779FF"), point.alpha = 0.7, background.indices = NULL,
    pch = 1, labels = c("sub-space", "total space"), legend = TRUE) {

    total_coors <- as.ppp(as.matrix(X[, dimensions]), c(range(X[, dimensions[1]]),
        range(X[, dimensions[2]])))
    total_space <- raster(density.ppp(total_coors))

    xlim <- range(X[, dimensions[1]])

    ylim <- if (is.null(background.indices))
        range(X[, dimensions[2]]) else c(min(X[, dimensions[2]]), c(max(X[, dimensions[2]]) + (max(X[, dimensions[2]]) -
        min(X[, dimensions[2]])) * 0.11))

    plot(x = X[, dimensions[1]], y = X[, dimensions[2]], col = "white", cex = basecex,
        xlab = dimensions[1], ylab = dimensions[2], xlim = xlim, ylim = ylim, cex.lab = basecex,
        xaxs = "i", yaxs = "i")

    # add background group
    if (!is.null(background.indices)) {

        # get points
        Y_bg <- X[background.indices, ]

        if (nrow(Y_bg) > 1 & point.alpha > 0)
            points(Y_bg[, dimensions[1]], Y_bg[, dimensions[2]], col = adjustcolor(colors[2],
                alpha.f = point.alpha), pch = pch)

        if (nrow(Y_bg) > 4) {
            coors <- as.ppp(as.matrix(Y_bg[, dimensions]), c(range(Y_bg[, dimensions[1]]),
                range(Y_bg[, dimensions[2]])))

            raster_dens <- raster(density.ppp(coors))

            palpre <- colorRampPalette(c("white", colors[2]))
            colsn <- 10
            palpre2 <- sapply(1:colsn, function(x) adjustcolor(col = palpre(colsn)[x],
                alpha = seq(0.1, 0.9, length.out = 20)[x]))

            image(raster_dens, add = TRUE, col = palpre2)
        }

    }

    Y <- X[indices, ]

    if (nrow(Y) > 1 & point.alpha > 0)
        points(Y[, dimensions[1]], Y[, dimensions[2]], col = adjustcolor(colors[1],
            alpha.f = point.alpha), pch = pch)

    if (nrow(Y) > 4) {
        coors <- as.ppp(as.matrix(Y[, dimensions]), c(range(Y[, dimensions[1]]),
            range(Y[, dimensions[2]])))

        raster_dens <- raster(density.ppp(coors))

        palpre <- colorRampPalette(c("white", colors[1]))
        colsn <- 10
        palpre2 <- sapply(1:colsn, function(x) adjustcolor(col = palpre(colsn)[x],
            alpha = seq(0.1, 0.9, length.out = 20)[x]))

        image(raster_dens, add = TRUE, col = palpre2)
    }

    usr <- par()$usr

    # legend
    if (!is.null(background.indices) & legend)
        legend(y = rep(usr[4] + ((usr[4] - usr[3]) * 0.05), 2), x = c(usr[1] + (usr[2] -
            usr[1]) * 0.25, usr[1] + (usr[2] - usr[1]) * 0.75), col = colors, legend = labels,
            pch = c(pch, pch), cex = basecex * 1.5, bty = "n", horiz = TRUE, xjust = 0)

    title(title, cex.main = basecex * 1.2)
}

# trait_space_plot(X = sels, indices = which(sels$ID == '395WBHM'), dimensions
# = c('TSNE1', 'TSNE2'), point.alpha = 0.2, pch = 1) trait_space_plot(X = sels,
# indices = which(sels$ID == '395WBHM'), dimensions = c('TSNE1', 'TSNE2'),
# background.indices = which(sels$ID != '395WBHM' & sels$record.group ==
# sels$record.group[sels$ID == '395WBHM'][1]), point.alpha = 0.2, pch = 1)
# Y <- sels[sels$ID == '395WBHM',] Y <- sels[sels$ID %in% c('395WBHM',
# '394WBHM', '360YGHM', '385YBHM'), ]

# trait_overlap(Y, dimensions = c('TSNE1', 'TSNE2'), group = 'ID', type =
# 'density', overlap.type = 'assymetric')


# X = Y dimensions = c('TSNE1', 'TSNE2') level <- '395WBHM' basecex <- 1 group
# <- 'ID' indices <- which(X[, group] == level)

trait_space_plot(X = sels, indices = which(sels$ID == "395WBHM"), dimensions = c("TSNE1",
    "TSNE2"), background.indices = which(sels$ID != "395WBHM" & sels$record.group ==
    sels$record.group[sels$ID == "395WBHM"][1]), point.alpha = 0.2, pch = 1)

# with points
out <- pblapply(unique(sels$ID), function(x) {

    tiff(filename = file.path("./images", paste(x, sels$treatment[sels$ID == x][1],
        "v2.tiff", sep = "-")), res = 120, width = 800, height = 800)

    par(mfrow = c(2, 2))

    # W <- sels[sels$ID == x, ]

    for (i in 1:4) try(trait_space_plot(X = sels, dimensions = c("TSNE1", "TSNE2"),
        indices = which(sels$ID == x & sels$week == i), background.indices = which(sels$ID !=
            x & sels$week == i & sels$record.group == sels$record.group[sels$ID ==
            x][1]), title = paste(x, "week", i), basecex = 1, labels = c(x, "group"),
        legend = FALSE), silent = TRUE)

    dev.off()
})

Session information

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] PhenotypeSpace_0.1.0  warbleR_1.1.26        NatureSounds_1.0.4   
##  [4] seewave_2.1.8         tuneR_1.3.3.1         vegan_2.5-7          
##  [7] lattice_0.20-44       permute_0.9-7         raster_3.5-15        
## [10] spatstat_2.3-3        spatstat.linnet_2.3-2 spatstat.core_2.4-0  
## [13] rpart_4.1-15          nlme_3.1-152          spatstat.random_2.1-0
## [16] spatstat.geom_2.3-2   spatstat.data_2.1-2   GGally_2.1.2         
## [19] adehabitatHR_0.4.19   adehabitatLT_0.3.25   CircStats_0.2-6      
## [22] boot_1.3-28           adehabitatMA_0.3.14   ade4_1.7-18          
## [25] deldir_1.0-6          sp_1.4-6              MASS_7.3-54          
## [28] Rtsne_0.15            randomForest_4.7-1    formatR_1.11         
## [31] knitr_1.37            kableExtra_1.3.4      ggplot2_3.3.5        
## [34] viridis_0.6.2         viridisLite_0.4.0     pbapply_1.5-0        
## [37] readxl_1.3.1          lubridate_1.8.0       remotes_2.4.2        
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_2.0-3      rjson_0.2.21          ellipsis_0.3.2       
##  [4] rstudioapi_0.13       proxy_0.4-26          farver_2.1.0         
##  [7] fansi_1.0.2           xml2_1.3.3            codetools_0.2-18     
## [10] splines_4.1.0         polyclip_1.10-0       jsonlite_1.8.0       
## [13] cluster_2.1.2         spatstat.sparse_2.1-0 compiler_4.1.0       
## [16] httr_1.4.2            assertthat_0.2.1      Matrix_1.3-4         
## [19] fastmap_1.1.0         cli_3.2.0             htmltools_0.5.2      
## [22] tools_4.1.0           gtable_0.3.0          glue_1.6.1           
## [25] dplyr_1.0.8           Rcpp_1.0.8            cellranger_1.1.0     
## [28] jquerylib_0.1.4       vctrs_0.3.8           svglite_2.1.0        
## [31] xfun_0.29             stringr_1.4.0         rvest_1.0.2          
## [34] lifecycle_1.0.1       goftest_1.2-3         terra_1.5-21         
## [37] scales_1.1.1          spatstat.utils_2.3-0  parallel_4.1.0       
## [40] RColorBrewer_1.1-2    yaml_2.3.5            gridExtra_2.3        
## [43] sass_0.4.0            reshape_0.8.8         stringi_1.7.6        
## [46] highr_0.9             rlang_1.0.1           pkgconfig_2.0.3      
## [49] systemfonts_1.0.4     dtw_1.22-3            bitops_1.0-7         
## [52] evaluate_0.15         purrr_0.3.4           tensor_1.5           
## [55] labeling_0.4.2        tidyselect_1.1.2      plyr_1.8.6           
## [58] magrittr_2.0.2        R6_2.5.1              fftw_1.0-6.1         
## [61] generics_0.1.2        DBI_1.1.1             pillar_1.7.0         
## [64] withr_2.4.3           mgcv_1.8-36           abind_1.4-5          
## [67] RCurl_1.98-1.6        tibble_3.1.6          crayon_1.5.0         
## [70] utf8_1.2.2            rmarkdown_2.11        grid_4.1.0           
## [73] digest_0.6.29         webshot_0.5.2         signal_0.7-7         
## [76] munsell_0.5.0         bslib_0.2.5.1